Trabajo Grupal 2

Jose Castañeda Chilon, Jose Joselito Carhuaricra Cusipuma, Joaquín Antonio Castañón Vilca

Librerías

Las librerías utilizadas son las siguientes:

library(olsrr)
library(BSDA)
library(MVN)
library(tidyverse)
library(hrbrthemes)
library(viridis)
library(ggplot2)
library(corrplot)
library(psych)
library(lmtest)
library(car)
library(tidymodels)
library(reshape2)
library(MLmetrics)
library(plotly)

Carga de Datos

data(trees)
attach(trees)
head(trees)
#>   Girth Height Volume
#> 1   8.3     70   10.3
#> 2   8.6     65   10.3
#> 3   8.8     63   10.2
#> 4  10.5     72   16.4
#> 5  10.7     81   18.8
#> 6  10.8     83   19.7

Preguntas

Pregunta 1

Construya gráficas apropiadas para las variables (histogramas, diagramas de dispersión por pares, boxplots para datos simétricos y asimétricos

par(mfrow=c(2,2))
boxplot(trees$Height, main='Height', col = 'red')
boxplot(trees$Girth, main='Girth', col = 'blue')
hist(trees$Height, main='Histogram of Height')
hist(trees$Girth, main='Histogram of Girth')
par(mfrow=c(1,1))

ggplot(trees, aes(x=Girth, y=Volume)) + geom_point(size=3) + geom_smooth(method = lm)
ggplot(trees, aes(x=Height, y=Volume)) + geom_point(size=3) + geom_smooth(method = lm)
ggplot(trees, aes(x=Height, y=Girth)) + geom_point(size=3)

cor(trees)
#>            Girth    Height    Volume
#> Girth  1.0000000 0.5192801 0.9671194
#> Height 0.5192801 1.0000000 0.5982497
#> Volume 0.9671194 0.5982497 1.0000000

pairs.panels(trees,
             smooth = TRUE,      # Si TRUE, dibuja ajuste suavizados de tipo loess
             scale = FALSE,      # Si TRUE, escala la fuente al grado de correlación
             density = TRUE,     # Si TRUE, añade histogramas y curvas de densidad
             ellipses = TRUE,    # Si TRUE, dibuja elipses
             method = "pearson", # Método de correlación (también "spearman" o "kendall")
             pch = 21,           # Símbolo pch
             lm = FALSE,         # Si TRUE, dibuja un ajuste lineal en lugar de un ajuste LOESS
             cor = TRUE,         # Si TRUE, agrega correlaciones
             jiggle = FALSE,     # Si TRUE, se añade ruido a los datos
             factor = 2,         # Nivel de ruido añadido a los datos
             hist.col = 4,       # Color de los histogramas
             stars = TRUE,       # Si TRUE, agrega el nivel de significación con estrellas
             ci = TRUE) 

Pregunta 2

Ejecute la prueba de bondad de ajuste a una normal trivariada. Interprete.

mvn(trees,mvnTest = c("mardia")) # NO
#> $multivariateNormality
#>              Test          Statistic            p value Result
#> 1 Mardia Skewness   20.9802564628517 0.0212316600572144     NO
#> 2 Mardia Kurtosis -0.163815893754234  0.869876079301746    YES
#> 3             MVN               <NA>               <NA>     NO
#> 
#> $univariateNormality
#>           Test  Variable Statistic   p value Normality
#> 1 Shapiro-Wilk   Girth      0.9412    0.0889    YES   
#> 2 Shapiro-Wilk  Height      0.9655    0.4034    YES   
#> 3 Shapiro-Wilk  Volume      0.8876    0.0036    NO    
#> 
#> $Descriptives
#>         n     Mean   Std.Dev Median  Min  Max  25th  75th       Skew   Kurtosis
#> Girth  31 13.24839  3.138139   12.9  8.3 20.6 11.05 15.25  0.5010559 -0.7109412
#> Height 31 76.00000  6.371813   76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846   24.2 10.2 77.0 19.40 37.30  1.0132739  0.2460393
mvn(trees,mvnTest = c("hz"))    # NO
#> $multivariateNormality
#>            Test      HZ    p value MVN
#> 1 Henze-Zirkler 0.92118 0.03216314  NO
#> 
#> $univariateNormality
#>           Test  Variable Statistic   p value Normality
#> 1 Shapiro-Wilk   Girth      0.9412    0.0889    YES   
#> 2 Shapiro-Wilk  Height      0.9655    0.4034    YES   
#> 3 Shapiro-Wilk  Volume      0.8876    0.0036    NO    
#> 
#> $Descriptives
#>         n     Mean   Std.Dev Median  Min  Max  25th  75th       Skew   Kurtosis
#> Girth  31 13.24839  3.138139   12.9  8.3 20.6 11.05 15.25  0.5010559 -0.7109412
#> Height 31 76.00000  6.371813   76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846   24.2 10.2 77.0 19.40 37.30  1.0132739  0.2460393
mvn(trees,mvnTest = c("royston")) # NO
#> $multivariateNormality
#>      Test        H   p value MVN
#> 1 Royston 8.331369 0.0170622  NO
#> 
#> $univariateNormality
#>           Test  Variable Statistic   p value Normality
#> 1 Shapiro-Wilk   Girth      0.9412    0.0889    YES   
#> 2 Shapiro-Wilk  Height      0.9655    0.4034    YES   
#> 3 Shapiro-Wilk  Volume      0.8876    0.0036    NO    
#> 
#> $Descriptives
#>         n     Mean   Std.Dev Median  Min  Max  25th  75th       Skew   Kurtosis
#> Girth  31 13.24839  3.138139   12.9  8.3 20.6 11.05 15.25  0.5010559 -0.7109412
#> Height 31 76.00000  6.371813   76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846   24.2 10.2 77.0 19.40 37.30  1.0132739  0.2460393
mvn(trees,mvnTest = c("dh")) # NO
#> $multivariateNormality
#>             Test        E df      p value MVN
#> 1 Doornik-Hansen 124.2377  6 2.096607e-24  NO
#> 
#> $univariateNormality
#>           Test  Variable Statistic   p value Normality
#> 1 Shapiro-Wilk   Girth      0.9412    0.0889    YES   
#> 2 Shapiro-Wilk  Height      0.9655    0.4034    YES   
#> 3 Shapiro-Wilk  Volume      0.8876    0.0036    NO    
#> 
#> $Descriptives
#>         n     Mean   Std.Dev Median  Min  Max  25th  75th       Skew   Kurtosis
#> Girth  31 13.24839  3.138139   12.9  8.3 20.6 11.05 15.25  0.5010559 -0.7109412
#> Height 31 76.00000  6.371813   76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846   24.2 10.2 77.0 19.40 37.30  1.0132739  0.2460393
mvn(trees,mvnTest = c("energy")) # NO
#> $multivariateNormality
#>          Test Statistic p value MVN
#> 1 E-statistic  1.191647   0.005  NO
#> 
#> $univariateNormality
#>           Test  Variable Statistic   p value Normality
#> 1 Shapiro-Wilk   Girth      0.9412    0.0889    YES   
#> 2 Shapiro-Wilk  Height      0.9655    0.4034    YES   
#> 3 Shapiro-Wilk  Volume      0.8876    0.0036    NO    
#> 
#> $Descriptives
#>         n     Mean   Std.Dev Median  Min  Max  25th  75th       Skew   Kurtosis
#> Girth  31 13.24839  3.138139   12.9  8.3 20.6 11.05 15.25  0.5010559 -0.7109412
#> Height 31 76.00000  6.371813   76.0 63.0 87.0 72.00 80.00 -0.3568773 -0.7233677
#> Volume 31 30.17097 16.437846   24.2 10.2 77.0 19.40 37.30  1.0132739  0.2460393

Como se puede observar, ninguna de las pruebas de bondad de ajuste mencionan que es los datos son Normales Multivariados, ya sea por Mardia, Henze-Zirkler, Royston, Doornik-Hansen y E-statistic. Sin embargo, por fines académicos se seguira trabajando con una regresión lineal múltiple.

Pregunta 3

Calcule las correlaciones con la función mixedcor de R

cor(trees)
#>            Girth    Height    Volume
#> Girth  1.0000000 0.5192801 0.9671194
#> Height 0.5192801 1.0000000 0.5982497
#> Volume 0.9671194 0.5982497 1.0000000
mixedCor(trees)
#> Call: mixedCor(data = trees)
#>        Girth Heght Volum
#> Girth  1.00             
#> Height 0.52  1.00       
#> Volume 0.97  0.60  1.00

Pregunta 4

Encuentre los estimadores puntuales y la ecuación de ajuste. Considere la variable “Volume” como la variable respuesta.

modelo <- lm(Volume~Girth + Height, data = trees)
summary(modelo)
#> 
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = trees)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.4065 -2.6493 -0.2876  2.2003  8.4847 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
#> Girth         4.7082     0.2643  17.816  < 2e-16 ***
#> Height        0.3393     0.1302   2.607   0.0145 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.882 on 28 degrees of freedom
#> Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
#> F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

\(Volume = -57.9877 + 4.7082*Girth + 0.3393*Height\)

Pregunta 5

Interprete el coeficiente de regresión asociada a la variable Girth

Por cada pulgada incrementada en la variable Girth, se reflejará como un aumento de 4.7082 pies cubicos en la variable respuesta “Volume”. La variable es significativa ya que tiene un valor p menor al valor de significancia de 0.05* (2e-16)

Pregunta 6

Presente la forma de calcular el Se. Interprete.

predicts <- cbind(predict(modelo, trees[,-3]))
originals <- c(trees$Volume)
se_table <- cbind(predicts, originals)
colnames(se_table) <- c('predicts', 'originals')
se_table_melt <- melt(se_table, id.vars='index', variable.name = 'series')
head(se_table)
#>    predicts originals
#> 1  4.837660      10.3
#> 2  4.553852      10.3
#> 3  4.816981      10.2
#> 4 15.874115      16.4
#> 5 19.869008      18.8
#> 6 21.018327      19.7
ggplot(se_table_melt, aes(Var1, value)) + geom_line(aes(colour = Var2))

Calculo del SE:

se_table <- data.frame(se_table)
se_table$SS <- ((se_table$originals-se_table$predicts)^2)
head(se_table)
#>    predicts originals         SS
#> 1  4.837660      10.3 29.8371621
#> 2  4.553852      10.3 33.0182211
#> 3  4.816981      10.2 28.9768907
#> 4 15.874115      16.4  0.2765548
#> 5 19.869008      18.8  1.1427790
#> 6 21.018327      19.7  1.7379860

SE <- sqrt(sum(se_table$SS)/(length(se_table$originals)-3))
print(SE)
#> [1] 3.881832
SE <- sqrt(sum(se_table$SS)/(length(se_table$originals)-3))
print(SE)
#> [1] 3.881832

Este valor también se puede obtener de la siguiente manera:

summary(modelo)
#> 
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = trees)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.4065 -2.6493 -0.2876  2.2003  8.4847 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
#> Girth         4.7082     0.2643  17.816  < 2e-16 ***
#> Height        0.3393     0.1302   2.607   0.0145 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.882 on 28 degrees of freedom
#> Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
#> F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

Residual standard error : \(3.882 = SE\)

Pregunta 7

Construya intervalos de confianza para el promedio del volumen y para un valor individual para Girth=11.5 y Height=73.

y <- mean(trees$Volume)

# promedio
tsum.test(mean.x=y, s.x=SE, n.x=length(se_table$originals), conf.level=0.95, alternative = 'two.sided')$conf
#> [1] 28.74710 31.59484
#> attr(,"conf.level")
#> [1] 0.95

# valor en especifico
new <- data.frame(cbind(11.5, 73))
predict(modelo, newdata = new, interval = 'prediction', level = 0.95)
#>          fit       lwr      upr
#> 1   4.837660 -3.561809 13.23713
#> 2   4.553852 -3.962908 13.07061
#> 3   4.816981 -3.809144 13.44311
#> 4  15.874115  7.690594 24.05764
#> 5  19.869008 11.451358 28.28666
#> 6  21.018327 12.469920 29.56673
#> 7  16.192688  7.797083 24.58829
#> 8  19.245949 11.092268 27.39963
#> 9  21.413021 13.103698 29.72234
#> 10 20.187581 12.047515 28.32765
#> 11 22.015402 13.775543 30.25526
#> 12 21.468465 13.327934 29.60900
#> 13 21.468465 13.327934 29.60900
#> 14 20.506154 12.270386 28.74192
#> 15 23.954110 15.854249 32.05397
#> 16 27.852203 19.760075 35.94433
#> 17 31.583966 23.126434 40.04150
#> 18 33.806482 25.303646 42.30932
#> 19 30.600978 22.388655 38.81330
#> 20 28.697035 19.945837 37.44823
#> 21 34.388184 26.295494 42.48087
#> 22 36.008319 27.878179 44.13846
#> 23 35.385260 27.237521 43.53300
#> 24 41.768998 33.386120 50.15188
#> 25 44.877702 36.655198 53.10021
#> 26 50.942868 42.647210 59.23853
#> 27 52.223751 43.899133 60.54837
#> 28 53.428513 45.064545 61.79248
#> 29 53.899329 45.522483 62.27618
#> 30 53.899329 45.522483 62.27618
#> 31 68.515305 59.707135 77.32347

Pregunta 8

Es la variable Height significativa? Ejecute la prueba de hipótesis completa.Presente la forma de calcular el p-valor.

H0 = No es significativo H1 = Es significtivo

anova(modelo)
#> Analysis of Variance Table
#> 
#> Response: Volume
#>           Df Sum Sq Mean Sq  F value  Pr(>F)    
#> Girth      1 7581.8  7581.8 503.1503 < 2e-16 ***
#> Height     1  102.4   102.4   6.7943 0.01449 *  
#> Residuals 28  421.9    15.1                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2*pt(2.607, 28, lower.tail = FALSE)
#> [1] 0.01447721

Dado los resultados, se puede decir que la variable Height es significativa, ya que el valor p es menor al valor de alpha (0.05) por esto no se rechaza la hipótesis nula

Pregunta 9

Es el modelo significativo? Ejecute la prueba de hipótesis completa.Presente la forma de calcular el p-valor.

H0 = No es significativo H1 = Es significtivo

g <- summary(modelo)
g$fstatistic
#>    value    numdf    dendf 
#> 254.9723   2.0000  28.0000
pf(254.9723,3,28, lower.tail = FALSE)
#> [1] 1.998157e-20

El modelo es significativo ya que el valor p calculado es menor al alpha establecido, por ende no se rechaza la hipótesis nula.

Pregunta 10

Presente la forma de calcular el R2 ajustado e interprete

El valor de R2 se puede obtener mediante el siguiente comando:

summary(modelo)
#> 
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = trees)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.4065 -2.6493 -0.2876  2.2003  8.4847 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
#> Girth         4.7082     0.2643  17.816  < 2e-16 ***
#> Height        0.3393     0.1302   2.607   0.0145 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.882 on 28 degrees of freedom
#> Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
#> F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

Pero también se puede calcular de la siguiente forma

r2 <- anova(modelo)
r2 <-1 - (r2$`Sum Sq`[3]/(sum(r2$`Sum Sq`)))
r2_adj <- 1 - (((1-r2)*(length(trees$Volume)-1))/(length(trees$Volume)-length(trees[,-3])-1))
print(r2)
#> [1] 0.94795
print(r2_adj)
#> [1] 0.9442322

Como resultado se tiene que el modelo tiene un 94.8% de que los datos se ajusten al modelo establecido con las 2 variables (Height y Girth) y así obtener el Volume

trees$predicted <- predict(modelo)
trees$residuals <- residuals(modelo)
ggplot(trees, aes(x = Height, y = Volume)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +     # regression line  
  geom_segment(aes(xend = Height, yend = predicts), alpha = .2) +      # draw line from point to line
  geom_point(aes(color = abs(residuals), size = abs(residuals))) +  # size of the points
  scale_color_continuous(low = "green", high = "red") +             # colour of the points mapped to residual size - green smaller, red larger
  guides(color = FALSE, size = FALSE) +                             # Size legend removed
  geom_point(aes(y = predicts), shape = 1) +
  theme_bw()

ggplot(trees, aes(x = Girth, y = Volume)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +     # regression line  
  geom_segment(aes(xend = Girth, yend = predicts), alpha = .2) +      # draw line from point to line
  geom_point(aes(color = abs(residuals), size = abs(residuals))) +  # size of the points
  scale_color_continuous(low = "green", high = "red") +             # colour of the points mapped to residual size - green smaller, red larger
  guides(color = FALSE, size = FALSE) +                             # Size legend removed
  geom_point(aes(y = predicts), shape = 1) +
  theme_bw()

Pregunta 11

Calcule los residuales. Interprete.

predicts <- cbind(predict(modelo, trees[,-3]))
originals <- c(trees$Volume)
se_table <- cbind(predicts, originals)
colnames(se_table) <- c('predicts', 'originals')
se_table <- data.frame(se_table)
se_table$residuals <- residuals(modelo)
head(se_table)
#>    predicts originals  residuals
#> 1  4.837660      10.3  5.4623403
#> 2  4.553852      10.3  5.7461484
#> 3  4.816981      10.2  5.3830187
#> 4 15.874115      16.4  0.5258848
#> 5 19.869008      18.8 -1.0690084
#> 6 21.018327      19.7 -1.3183270

El como se puede observar, existe un grado de cuanto a la predicción del Volume con el modelo, principalmente con los 3 primeros datos, pudiendo ser estos outliers. Para obtener un mejor modelo se recomienda poder extraer una mayor cantidad de datos.

H0 = Los residuales son homocedasticos, varianza constante H1 = Los residuales no homocedasticos (Heterosedasticos)

Graficar los residuos frente a los valores ajustados por el modelo, e identificar si existe un patrón cónico u otro patrón. Idealmente deberían distribuirse de forma aleatoria en torno a 0. También podemos recurrir al test de Breusch-Pagan como contraste de homocedasaticidad.

#no estandarizada
plot(modelo, which=1, col=c("red"))

bptest(modelo)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  modelo
#> BP = 2.4681, df = 2, p-value = 0.2911

valor p = 0.2911 > 0.05, no se rechaza la H0, es decir que los residuos tienen varianza constante

Pregunta 12

Calcule los residuales estandarizados y estudentizados. Interprete.

#estandarizada y estudientizado
se_table$residuals_standard <- rstandard(modelo)
se_table$residuals_studentized <- rstudent(modelo)
# standarized graph
plot(modelo, which=2, col=c("red"))
ols_plot_resid_stand(modelo)

# studentized graph
hist(se_table$residuals_studentized, freq = FALSE)
qqPlot(modelo)
#> [1] 18 31
ols_plot_resid_stud(modelo)

Lo que puede observar es que los valores residuales estandarizados y estudiantizados tienden a acercarse a la linea de 45º, por lo que, podemos decir que los valores de los residuales poseen una distribución normal

Pregunta 13

Presente la forma de calcular la distancia de Cook de los cinco primeros registros. Interprete.

ols_plot_cooksd_bar(modelo)

ols_plot_resid_lev(modelo)

Se puede observar que los registros 20 y 31 son outliers

cooksD <- cooks.distance(modelo)
se_table$cooksD <- cbind(cooksD)
head(se_table)
#>    predicts originals  residuals residuals_standard residuals_studentized
#> 1  4.837660      10.3  5.4623403          1.4964901             1.5320694
#> 2  4.553852      10.3  5.7461484          1.6029462             1.6516683
#> 3  4.816981      10.2  5.3830187          1.5284555             1.5677398
#> 4 15.874115      16.4  0.5258848          0.1396700             0.1372010
#> 5 19.869008      18.8 -1.0690084         -0.2936751            -0.2888284
#> 6 21.018327      19.7 -1.3183270         -0.3696163            -0.3638447
#>         cooksD
#> 1 0.0977927647
#> 2 0.1478462779
#> 3 0.1673192079
#> 4 0.0004091116
#> 5 0.0039449244
#> 6 0.0084012068
influential <- cooksD[(cooksD > (3*mean(cooksD, na.rm = TRUE)))]
influential
#>         3        18        31 
#> 0.1673192 0.1775359 0.6052326
cooks_table <- se_table

Calculo manual:

Hmat <- hatvalues(modelo)
cooks_table$h <- c(Hmat)
data(trees)
attach(trees)
head(trees)
#>   Girth Height Volume
#> 1   8.3     70   10.3
#> 2   8.6     65   10.3
#> 3   8.8     63   10.2
#> 4  10.5     72   16.4
#> 5  10.7     81   18.8
#> 6  10.8     83   19.7
k <- length(trees[,-3])
cooks_table$cookD_calc <- ((cooks_table$residuals_standard^2)*cooks_table$h)/((k+1)*(1-cooks_table$h)) # k+1 => k: numero de variables predictoras

En la tabla cooks_table, en la columna “cookD_calc” se muestran los resultados de la distancia de cook que antes se calculo con la formula:

head(cooks_table, 5)
#>    predicts originals  residuals residuals_standard residuals_studentized
#> 1  4.837660      10.3  5.4623403          1.4964901             1.5320694
#> 2  4.553852      10.3  5.7461484          1.6029462             1.6516683
#> 3  4.816981      10.2  5.3830187          1.5284555             1.5677398
#> 4 15.874115      16.4  0.5258848          0.1396700             0.1372010
#> 5 19.869008      18.8 -1.0690084         -0.2936751            -0.2888284
#>         cooksD          h   cookD_calc
#> 1 0.0977927647 0.11582883 0.0977927647
#> 2 0.1478462779 0.14720958 0.1478462779
#> 3 0.1673192079 0.17686186 0.1673192079
#> 4 0.0004091116 0.05919131 0.0004091116
#> 5 0.0039449244 0.12066468 0.0039449244
#Pregunta 14

Presente detalladamente la forma de calcular la matriz variancias-covariancias del vector de estimadores de beta.

matrix_trees <- trees[, 1:2]
matrix_trees$r <- c(rep(1, 31))
matrix_trees <- matrix_trees[, c(3,1,2)]
matrix_trees <- as.matrix(matrix_trees)
vcov(modelo)
#>             (Intercept)       Girth      Height
#> (Intercept)  74.6189461  0.43217138 -1.05076889
#> Girth         0.4321714  0.06983578 -0.01786030
#> Height       -1.0507689 -0.01786030  0.01693933
(SE^2)*solve(t(matrix_trees)%*%(matrix_trees))
#>                 r       Girth      Height
#> r      74.6189461  0.43217138 -1.05076889
#> Girth   0.4321714  0.06983578 -0.01786030
#> Height -1.0507689 -0.01786030  0.01693933

Pregunta 15

Calcule el VIF. Comente sus resultados

vif(modelo)
#>   Girth  Height 
#> 1.36921 1.36921

El factor de inflación de varianza (VIF) muestra que tanto para la variable Girth como Height toma valores menores a 10, lo que significa que ambas variables entran al modelo.

Pregunta 16

Elija al azar 70% de los datos y proceda a ejecutar la RLM para obtener la ecuación de ajuste.

division <- trees %>% initial_split(prop = 0.7, strata = Volume)
training_set <- division  %>%  training()
testing_set <- division  %>%  testing() 

length(training_set$Volume)
#> [1] 21
length(testing_set$Volume)
#> [1] 10

Desarrollo de RLM:

modelo_train <- lm(Volume ~ Girth + Height, data = training_set)
summary(modelo_train)
#> 
#> Call:
#> lm(formula = Volume ~ Girth + Height, data = training_set)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.1422 -2.1230 -0.3601  2.9643  5.1460 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -48.9374    13.2267  -3.700  0.00164 ** 
#> Girth         4.6903     0.3118  15.043 1.23e-11 ***
#> Height        0.2267     0.1885   1.203  0.24469    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.693 on 18 degrees of freedom
#> Multiple R-squared:  0.9447, Adjusted R-squared:  0.9385 
#> F-statistic: 153.6 on 2 and 18 DF,  p-value: 4.864e-12

Pregunta 17

Con la ecuación obtenida en el paso anterior, calcule los valores ajustados del 30% restante

predicts <- cbind(predict(modelo, testing_set[, 1:2]))
predicts_table <- cbind(predicts, c(testing_set[,3]))
colnames(predicts_table) <- c('predicts', 'originals')
predicts_table <- as.data.frame(predicts_table)
predicts_table
#>     predicts originals
#> 3   4.816981      10.2
#> 4  15.874115      16.4
#> 6  21.018327      19.7
#> 7  16.192688      15.6
#> 14 20.506154      21.3
#> 19 30.600978      25.7
#> 20 28.697035      24.9
#> 22 36.008319      31.7
#> 30 53.899329      51.0
#> 31 68.515305      77.0

Pregunta 18

Calcule el MSE del 30% restante.

mean((predicts_table$original-predicts_table$predicts)^2)
#> [1] 16.93677
MSE(predicts_table$predicts, predicts_table$original)
#> [1] 16.93677

Pregunta 19

Explique la diferencia entre el cuadrado medio del error del análisis de variancia y el MSE

Según los resultados del anova el MSE es 12.1 (es el resultado de dividir la suma de cuadrados del error entre n - k - 1) donde (k es numero de coeficientes = 3, n es el tamaño de la muestra = 21), mientras que el MSE calculado resultó 13.47 (es el resultado de la división entre n).

Pregunta 20

Grafique los datos en 3D usando R

fig <- plot_ly(x = trees$Girth, y = trees$Height, z = trees$Volume)
fig